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ABSTRACT 

We present a new approach to foreground removal for Cosmic Microwave Back- 
ground (CMB) maps. Rather than relying on prior knowledge about the foreground 
components, we first extract the necessary information about them directly from the 
microwave sky maps by taking differences of temperature maps at different frequencies. 
These difference maps, which we refer to as internal templates, consist only of linear 
combinations of galactic foregrounds and noise, with no CMB component. We obtain 
the foreground cleaned maps by fitting these internal templates to, and subsequently 
subtracting the appropriately scaled contributions of them from, the CMB dominated 
channels. The fitting operation is performed in wavelet space, making the analysis fea- 
sible at high resolution with only a minor loss of precision. Applying this procedure 
to the WMAP data, we obtain a power spectrum that matches the spectrum obtained 
by the WMAP team at the signal dominated scales. The fact that we obtain basi- 
cally identical results without using any external templates has considerable relevance 
for future observations of the CMB polarization, where very little is known about the 
galactic foregrounds. Finally, we have revisited previous claims about a north-south 
power asymmetry on large angular scales, and confirm that these remain unchanged 
with this completely different approach to foreground separation. This also holds when 
fitting the foreground contribution independently to the northern and southern hemi- 
sphere indicating that the asymmetry is unlikely to have its origin in different foreground 
properties of the hemispheres. This conclusion is further strengthened by the lack of 
any observed frequency dependence. 

Subject headings: (cosmology:) cosmic microwave background — cosmology: observa- 
tions — methods: data analysis — methods: statistical 

1. Introduction 

The Wilkinson Microwave Anisotropy Probe ( WMAP) satellite (Bennett et al. 2003a) and other 
recent ground-based and balloon-borne experiments have provided high sensitivity observations of 

the Cosmic Microwave BackgroTind (CMB). High-precision estimates of the cosmological parameters 
have been derived from these data, and these will improve yet further with data from Planck and 
near-future ground based experiments. Later still, high sensitivity polarization data of the CMB 
will be available and the physics of the early universe can be studied with even higher fidelity. 
However, in order to obtain reliable estimates of the cosmological parameters, control of systematic 
effects and the different sources of foreground contamination is of the utmost importance. In this 
paper we address the latter problem. 

There are three well-understood sources of diffuse galactic contamination. Synchrotron emis- 
sion from the Galaxy originates in relativistic cosmic ray (CR) electrons spiraling in the Galactic 
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magnetic field. Free-free emission is the bremsstrahlung radiation resulting from the Coulomb in- 
teraction between free electrons and ions in the Galaxy. Thermal dust emission arises from grains 
large enough to be in thermal equilibrium with the interstellar radiation field. Evidence of an 
additional component initially arose from the cross-correlation of the COBE-DMR data with the 
DIRBE map of thermal dust emission at 140//m (Kogut et al. 1996) that revealed anomalous emis- 
sion with a spectrum rising as the frequency decreased from 53 to 31.5 GHz. Both Banday et al. 
(2003) and Bennett et al. (2003b) suggested that this component was characterized by a power-law 
frequency spectrum with index -2.5 over the range ^^20 - 60 GHz. However, it remains unclear 
what the physical origin of such emission is - Bennett et al. (2003b) have proposed that it arises 
from hard synchrotron emission in star-forming regions close to the Galactic plane, whilst Draine 
Sz Lazarian (1998) suggest rotational emission from very small grains or 'spinning dust'. Recently, 
Watson et al. (2005) have shown that observations of the Perseus molecular cloud between 11 and 
17 GHz and augmented with the WMAP data can be adequately fitted by a spinning dust model. 
Nevertheless, we will refer to the putative component as anomalous dust. 

Fortunately, the foreground components have frequency spectra very different from the CMB, 
although the exact shapes of these spectra are not known. Most data on galactic emission have been 
taken at frequencies others than those used for CMB observations and it is not clear whether using 
these as templates for CMB foreground subtraction is a valid approach. Certainly spatial variations 
in the frequency dependence of the foregrounds will cause the templates to increasingly diverge 
from the true structure at microwave wavelengths. Nevertheless, given the lack of other reliable 
approaches, a template fitting and subtraction procedure was applied on the WMAP data (Bennett 
et al. 2003a) using external templates (Bennett et al. 2003b; Hinshaw et al. 2003). The Spectral 
Matching Independent Component Analysis (SMICA) (Delabrouille et al. 2003) method has been 
applied to the WMAP maps after they had already been cleaned using the standard external 
templates and evidence for a residual galactic component was found (Patanchon et al. 2005). The 
WMAP data were also corrected using the Internal Linear Combination (ILC) methods (Bennett 
et al. 2003b; Tegmark et al. 2003; Eriksen et al. 2004a) but unfortunately these methods are 
biased (Eriksen et al. 2004a). Other methods that have been developed for foreground subtraction 
(Barreiro et al. 2004; Hobson et al. 1998; Stolyarov et al. 2002, 2005; Maino et al. 2002, 2003; 
Donzelh et al. 2005; Brandt et al. 1994; Eriksen et al. 2005) have yet to be applied to the WMAP 
data. 

The importance of being able to make reliable foreground corrections on large angular scales 
with limited knowledge about the foreground components becomes even more apparent when con- 
sidering future CMB polarization observations. A very important test of inflation depends on high 

precision measurements of the predicted B-mode polarization of the CMB on large angular scales. 
However, very little is known about the polarization of the galactic emission. Future high sensitivity 
measurements of the CMB polarization arc likely to be highly dependent on foreground subtraction 
methods that do not make strong assumptions about the nature of the galactic emission. 

Here we propose a new foreground subtraction method for which no knowledge about the 
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morphology or frequency spectra of galactic components are required. The method does not require 
constant (in frequency) spectral indexes for the galactic emission components and will therefore 
not be biased by uncertainties in the changes of the spectral indexes over large frequency ranges. It 
will however require the spectral indexes to be constant in space over a given patch on the sphere. 
In the implementation presented here, we assume the spectral indexes to be constant over the full 
sphere or over hemispheres but the extension to smaller patches will be discussed. As the method 
increases the noise level in the data, it is at present mainly useful for studies of the larger scales 
where noise is not dominant. This is particularly relevant given that the largest angular scales 
in the WMAP data have indicated the presence of several anomalies, in particular a north-south 
asymmetry in the power spectrum (Eriksen et al. 2004b; Hansen et al. 2004). This result will be 
re-examined here with our new foreground corrected map. 

The method presented here takes advantage of the fact that the observed data already pro- 
vide information about the foregrounds. In particular, by computing differences of maps observed 
at different frequencies, (noisy) linear combinations of the foregrounds components are obtained. 
Fitting and subtracting such templates that are simple linear combinations of the foregrounds is 
mathematically equivalent to fitting and subtracting templates of the physical foreground compo- 
nents. The advantage of such linear combinations is that they can be obtained directly from the 
microwave observations themselves and therefore do not rely on other experiments. As a conse- 
quence, the foreground morphologies are likely to be well traced even in the presence of modest 
departures from a single spectral index in a given region. 

In §2, we present basis of the method and demonstrate its application in detail. Then, in §3 
and §4 the results of this procedure as applied to both simulated inputs and the WMAP sky maps 
are presented. Finally, in §5, we discuss how the method can be improved in the future. 



2. Methodological Basis: Fitting and Subtracting Internal Templates 

A pixelized temperature map of microwave observations can be written (in thermodynamic 
temperatures) as 

Nt 

Tr = Tf^ + n^ + Y,c^4- (1) 
t=i 

Here T,!^ is the observed temperature in pixel i for frequency channel u, J^cmb ^Yiq frequency 
independent CMB component, nj' is the instrumental noise and finally s* is the contribution from 
galactic foreground component t for pixel i. We assume a total of Nt different foreground compo- 
nents. The coefficients give the amplitude of the given foreground component t in channel u. In 
this paper we assume spatially constant foreground frequency spectra, as was also assumed by the 
WMAP team in generating the publicly available foreground cleaned maps. However, as discussed 
in the conclusions, the method can be extended to take into account variations of spectral indexes 
across the sky, but implementation of this is deferred to a future publication. 
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2.1. Internal templates 

The WMAP team used external templates obtained from observations at other frequencies as 
tracers of the foreground components s\. The corresponding coefficients were then obtained by 
a fitting procedure (the details of which are not clear) and the inferred foreground contributions 
subsequently sTibtracted from the sky maps. Here, we will not use external templates, but instead 
take the difference between two frequency channels. In such difference maps, which we call internal 
templates, the CMB cancels out and the remainder is a linear combination of foreground components 
plus instrumental noise. 

We can write the internal templates D'(^' , as 

Nt 

Df n' = ^(cr - cf )4 + 5nf, (2) 

t=i 

where 5nf^'^' = — nf is the noise of the internal template. Note that we have assumed the maps 
at all frequencies to have been smoothed to the same resolution. 

The natural physical 'basis' for the foreground components may be identified with dust, syn- 
chrotron and free-free emission. The internal templates are transformations into a new basis con- 
sisting of linear combinations of the physical components. Fitting and subtracting foreground 
components in this new basis is clearly mathematically equivalent to using the 'natural basis', pro- 
vided that the number of internal templates is equal to the number of components. Therefore, 
using the internal templates — (5nJ"^' as our new basis, equation (1) can be written as 

Nt 

Tr = tP'^'' + < + E ~^tiDi - Sn'), (3) 
t=i 

where the sum is performed over internal templates t = uv' for Nt different combinations of chan- 
nels. 

It is important to note that because the internal templates are noisy, subtracting them from 
the data will increase the overall noise level of the final map and significantly increase the error 
bars for power spectrum estimates at the highest multipoles. 



2.2. Fitting procedure 

When using external templates s\ for the different foreground components, the fitting procedure 
is performed by minimizing the x^, 

= E E(^'^ - E <S^CrjiT^ - E <4) (4) 

ij V t t 
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for the coefficients c^. Here tfie sum is performed over pixels i and j as well as channels i' including 
only the channels we wish to clean. The correlation matrix is defined by Cij = (TiTj) — {Ti){Tj) 
which can be obtained by Monte-Carlo simulations assuming a reasonable CMB spectrum. In this 
paper we use a similar fitting procedure, replacing sj — > Dj. In principle we should have used 
Dj — 8n\ (see equation 3) , but since the noise is known only in a statistical sense we are forced to 
use D\ and correct for the resulting noise bias. This bias correction procedure will be described in 
detail in Appendix A. 

In practise, the evaluation of equation (4) is not feasible for high resolutions. The covariance 
matrix is too large to be stored in a computer and too time consuming to invert. For this reason 
we adopt a wavelet space implementation of the above procedure. The advantage is that we 
can neglect pixel-pixel correlations, taking only into account the scale-scale correlations in the 
correlation matrix without significant loss of precision. 

Wc adopt the wavelet fitting procedure of, e.g., Vielva et al. (2004a) and Hansen et al. (2005) 
and define the cross-correlation coefficients 

Xf^Y.<s^ls, X'^' ^Y.'^ls^'s- (5) 

i i 

Here w'(g is the wavelet transform of for pixel i, channel u and scale 5, and w\g is the wavelet 
transform of the internal template D\. We use Spherical Mexican Hat wavelets (Martinez-Gonzalez 
et al. 2002) which are suitable for both template fitting (Vielva et al. 2004a; Hansen et al. 2005), 
point source detection (Vielva et al. 2003) and tests of non-Gaussianity (Vielva et al. 2004b; Mukher- 
jee & Wang 2004; Cabella et al. 2004). 

In wavelet space, the to minimize has the form 

= E E(^5^ - E -tXf)c-sUx^sf - E -T^F), (6) 

uT SS' t t 

where the summation runs over all target channels u, all internal templates (foreground compo- 
nents) T and all wavelet scales S and S' . The scale-scale correlation matrix Css' is obtained by 
Monte-Carlo simulations. Minimizing the x^, the best fit are given by a set of equations 



Y,SiMtf = B) (7) 

t 

for each template / and each frequency v. Here the matrix Mtf is given by 

^*/ = EE4^^5i'^F- (8) 

T SS' 

and the vector B'^ by 

Bf = J2J2^l^Cs's,X^ (9) 

T SS' 
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An estimate of the error bars on the coefficients can be obtained performing Monte-Carlo 
simulations of CMB, noise and 'foregrounds'. In these simulations, the internal templates obtained 
from the data arc used as 'foregrounds' with the best estimate coefficients obtained from the data. 
An additional correction to the Xg is necessary in this case as described in Appendix A. The error 
bars obtained in this manner has been shown to agree well with the real error bars. 

The wavelet scales to use in the analysis are chosen such that the scale-scale correlation matrix 
is well conditioned. In order to determine these scales, we adopt the following procedure: 

• Use MC simulations to obtain Cgs' for a large set of scales. 

• Define a limit a, say, 0.95. Start by the smallest allowed scale S and find the next scale by 
identifying at which scale the normalized covariance matrix Cgs' I V CssCs'S' has fallen to a. 

• Repeat the above procedure until the largest scale has been reached. Check whether the final 
correlation matrix is well conditioned. If not, decrease a and repeat. 

As in pixel space, a correction procedure for the noise bias has to be implemented. In Appendix 
A we will describe this procedure in detail as well as a procedure for estimating the level of the 
remaining bias. 

3. Testing the Method on Simulations 

The foreground subtraction procedure described above, which we will call Wavelet based high 
resolution Fitting of Internal Templates (WI-FIT), has been extensively tested on simulated maps. 
We generated a set of 500 Monte-Carlo simulations of CMB and noise, using the best-fit WMAP 
power-law power spectrum and noise properties corresponding to the five WMAP channels^. To 
these simulations, we added contributions from known galactic foregrounds based on specific tem- 
plates. In particular, for thermal dust we use the template provided by Schlcgcl ct al. (1998) 
extrapolated in frequency by Finkbeincr ct al. (1999), for synchrotron we adopt the 408 MHz sky 
survey of Haslam et al. (1982), and for the free-free contribution we assume that the emission can 
be traced by a template of emission (Finkbeiner 2004). The templates are then scaled to the 
WMAP frequencies using the weights given in Table 4 of Bennett et al. (2003b). These dust weights 
effectively include an anomalous dust contribution assuming that the putative emission can also be 
well traced by the thermal dust template. 

In order to validate the wavelet basis of our high resolution analysis, we first make an explicit 
comparison of the method to a full pixel space foreground subtraction procedure at lower resolution 
where it remains feasible. The maps were therefore smoothed to a common resolution of 5.5° and 



^These may be obtained from the Lambda website http://lambda.gsfc.nasa.gov/. 
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degraded to HEALPix^ resolution Nsibe = 32. Unless otherwise stated, the WMAP Kp2 sky cut 
was applied in all the analyses presented in this paper. Using the procedure with internal templates 
as described above, then after including the bias correction wc obtained unbiased estimates of the 
foreground coefficients in pixel space. Repeating the same procedure in wavelet space yields 
approximately 2 — 3% larger error bars showing that the loss in precision using this approach 
instead of the full pixel space procedure is small. Therefore, since analysis at higher resolution is 
in any case unfeasible in pixel space, we adopt the wavelet approach in what follows. 

We then attempted to calibrate the efficiency of internal versus external template fitting. The 
WI-FIT procedure was applied to maps of 1 degree FWHM resolution at A^side = 256 and compared 
to wavelet-based fitting of external templates. The results derived using 13 wavelet scales is shown 
in Fig. 1. The two upper maps show the Icr residuals for the external and internal template 
subtraction method when the internal templates were noise-free. These maps were obtained by 



where the sum is performed over all simulations s, in total Ng = 500, M? is the cleaned map for 
pixel i and realization s and T? is the input CMB realization. 

Unfortunately, the exercise is not particularly revealing, and mostly emphasises that external 
templates are likely to outperform internal templates when we know that they accurately describe 
the foregrounds present in the data. In the example shown here, the internal templates do better 
close to the galactic plane, but worse outside. Additionally, since the internal templates will 
be noisy, the la errors will be further infiated as shown in the lower map. Furthermore, the 
exact level of residuals will in general depend on the choice of internal templates (the specific 
difference combinations adopted in the analysis) as well as the nature of the foregrounds. In fact, 
the comparison somewhat violates the philosophy underpinning the WI-FIT method. The external 
templates arc measured at frequencies different from those used for studies of the CMB and are 
therefore not likely to accurately trace the foregrounds in the CMB dominated channels (although 
they do by construction in this test). The internal templates however, are linear combinations of 
the foregrounds at the channels used for CMB studies. The comparison is therefore not completely 
fair and the results shown in the plot do not adequately reflect a true comparison. We nevertheless 
present the results for completeness. 



^http://healpix.jpl. nasa.gov/ 
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Fig. 1. — Results from 500 Monte-Carlo simulations: The upper plot shows the la residuals y {5f) 
on the V-band from the wavelet fitting procedure using external templates. The middle and lower 
plot show the same result when using internal templates without and with noise respectively. Note 
that in these simulations the external templates are equal to the foregrounds added to the simulated 
maps. In reality, templates taken at frequencies different from the ones dominated by the CMB will 
differ from the foregrounds contaminating the CMB. This additional error is not reflected in the 
upper plot. Note further that the last case generally has higher noise variance than the two other 
cases, but the noise variance has not been included in the plots in order to make the foreground 
residuals clearer. 
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4. Application to the WMAP Data 
4.1. The maps 

Wc have applied the WI-FIT method to the first-year WMAP data smoothed to 1 degree 
FWHM at A'^siDE = 256 and using the Kp2 pixel mask. For each of the WMAP sky maps covering 
five frequency bands - K, Ka, Q, V and W (when multiple maps are available at the same frequency 
we take simple averages of the maps after convolution to the 1 degree FWHM beam) - three 
internal templates were constructed using the four remaining bands. For example, for the Q-band 
the internal templates K-Ka, Ka-V and V-W were generated, fitted to the Q-band sky map, and 
subtracted according to the above prescription. 

In Fig. 2, we show the cleaned K, Ka, Q, V and W maps. The template coefficients for the 
K band are so large that the noise level of the resultant cleaned map is too high to be of practical 
value in any further analysis and will not be considered in what follows. The maps shown can be 
obtained by combining the different WMAP bands according to the weights given in table 1. In 
the upper plot of Fig. 9 we also show a noise-weighted combination of the Q, V and W bands. We 
do not present error bars here as these error bars only are valid when the model (that the spectral 
indexes of the foregrounds are the same in all directions) is valid. We know that the assumption 
of spatially constant spectral indexes is wrong and the error bars would hence be misleading. 

In Fig. 2 there are no obvious differences between the Ka, Q, V and W bands. In order to 
check for possible foreground residuals, we have taken the differences between these sky maps. Since 

these foreground residuals are generally smaller than the instrumental noise level, the difference 
maps resemble pure noise without further processing. For visualisation purposes, we have applied 
a median filter with a 3 degree radius to suppress the noise on pixel scales. Fig. 3 shows the 
the difference maps for Q-V and V-W. For a perfect foreground subtraction, the difference maps 
should show little coherent structure beyond that expected for median-filtered pure noise. It is 
apparent that some small residuals remain outside the Kp2 cut. For comparison, the corresponding 
differences of the WMAP maps cleaned by External Template Fitting (ETF) provided by the 
WMAP collaboration are also shown. The ETF maps seem to have stronger residuals than the 
WI-FIT maps close to the galactic cut. The WI-FIT maps show stronger fluctuations over the 
whole sky, but this can partly be explained by the higher noise level in these maps. 

In Figure (4) we show the noise-corrected spectra of these four difference maps. For the Q-V 
difference, WI-FIT shows stronger residuals in the lowest multipole bin, but ETF has residuals well 
above £ = AO where the WI-FIT difference spectrum is consistent with pure noise. For the V-W 
map however, ETF shows stronger residuals in the smallest multipole bin whereas for all other bins 
the two maps are consistent. 

Note further that the Q-V ETF difference map shows strong similarities to the residual com- 
ponent obtained by the SMICA method. Patanchon et al. (2005) found evidence for a residual 




Fig. 2. — WI-FIT cleaned WMAP maps: The 5 WMAP channels smoothed to a common resolution 
of 1 degree FWHM and cleaned with the WI-FIT procedure applied outside the Kp2 galactic cut. 




Fig. 3. — The plots in the left column show the median filtered (3 degree radius top-hat) differences 
Q-V and V-W for the WI-FIT cleaned WMAP maps where only that part of the sky outside the Kp2 
mask was used in the fitting procedure. The plots in the right column show the same differences 
for the maps that have been foreground corrected by the WMAP collaboration using external 
templates. 
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Fig. 4. — The binned power spectra of the difference maps Q-V and V-W for the ETF and WI-FIT 
maps (The Kp2 galactic cut was used in the analysis). Black points: WI-FIT Q-V, green points: 
ETF Q-V, red points: WI-FIT V-W, blue points: ETF V-W. The error bars are obtained from 
pure noise simulations. 
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galactic component (their Figure 5) mainly present in the Q band but not detected in the V band^. 
It is therefore not unexpected that this component shows up when taking the difference Q-V. 

In the V-W difference maps one can in both cases see a blue region following the border of the 
galactic cut (this is more prominent in the ETF maps). As the W band is mainly contaminated by 
thermal dust, this can be interpreted as a sign of residual thermal dust contamination. That there 
is a residual thermal dust contamination in the W band is supported by table 1 where we have 
assumed a spectral index for each foreground component {(3s = —2.7 for synchrotron, (3^ = —2.15 
for free-free, (3^ = 2.2 for thermal dust and /?sd = —2.5 for anomalous dust) and estimated the 
proportions of the residual foreground components (see Eriksen et al. (2004a) for details). These 
numbers indicate the fraction of residual foreground contribution for each type relative to the 
uncorrected level at an adopted reference frequency (22.8 GHz for synchrotron and anomalous 
dust, 33.0 GHz for free- free and 93.5 GHz for thermal dust). 

In Fig. 5 wc show the filtered difference maps between the WTFIT and ETF maps and the 
ILC map provided by the WMAP collaboration. The Figures suggest the presence of rcsidiials 
at the 10% level compared to the CMB over the full sky in either the WI-FIT maps, ETF maps 
or both. Note that the WMAP-F,TF difference for the Q band shows similarities to the residual 
component detected with SMICA. The fact that the hot and cold areas are interchanged indicates 
that the residual is not present in the WI-FIT cleaned maps. A ring-like structure on large angular 
scales is also in all channels. It is not clear whether this residual structure can be attributed to 
the ETF or the WI-FIT cleaned maps. As the residual is positive in WI-FIT minus ETF, it can 
either be a strucure which is corrected for in ETF but not in WI-FIT, or an overcorrection in ETF. 
Looking at the power spectra of these difference maps (Fig. 6), the residual does not significantly 
affect the estimated CMB power spectra. 



4.2. The power spectrum 

In this section, we consider the implications of these residuals for power spectrum estimation. 
The MASTER algorithm (Hivon et al. 2002) was applied to the WI-FIT foreground cleaned maps 
outside the Kp2 sky cut. Before obtaining the spectra, the best fit mono- and dipole were subtracted 
from all maps. Fig. 6 presents the results (binned version in Fig. 7). For comparison, the power 
spectrum estimated by the WMAP collaboration (Hinshaw et al. 2003) from the ETF cleaned maps 
using the cross power spectrum is shown by a solid black line. The WI-FIT spectra are shown as 
green (Q), blue (V) and red (W) lines and the Ka spectrum is delineated by crosses. 

There is good agreement between the V and W band spectra estimated with the WI-FIT blind 
approach and the ETF cleaning method based on external observations of the galaxy at very dif- 
ferent frequencies. Small differences in the Q band spectrum probably indicate residual foreground 



^This feature is actually faintly visible in Figure 11 of Bennett et al. (2003b). 
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Fig. 5. — The plot shows the median filtered (3 degree radius top-hat) differences between the 
WI-FIT and ETF/ILC maps obtained by the WMAP collaboration for the three frequencies Q, V 
and W. 
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Fig. 6. — The power spectra of the WI-FIT cleaned WMAP maps for different bands compared to 
the power spectrum obtained by the WMAP team (soHd black line). Green line shows the result 
for the Q band, blue line for the V band, red line for the W band and black crosses for the Ka band. 
The lower green and red lines show the power spectra of the difference WI-FIT minus ETF for the 
Q and W band respectively. The power spectrum of the difference map Q-W for the WI-FIT maps 
is shown as a black line. 
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Fig. 7. — The binned power spectra of the WI-FIT cleaned WMAP maps for different bands 
compared to the power spectrum obtained by the WMAP team (histogram and error bars). Green 
crosses show the result for the Q band, blue crosses for the V band, red crosses for the W band 
and pink crosses for the Ka band. 
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contamination, more obviously present in the results for the highly foreground contaminated Ka 
channel. Nevertheless, both agree reasonably well with the V and W band spectra. Finally, we also 
show noise-corrected spectra of the WI-FIT and ETF difference maps. The foreground residuals 
outside the Kp2 cut which show up when taking the difference between the WI-FIT cleaned Q, V 
and W maps (Fig. 3) as well as the difference between the ETF and WI-FIT maps (Fig. 5) are 
clearly unlikely to cause significant differences in the CMB power spectrum. 

4.3. Pow^er asymmetries 

Eriksen et al. (2004b) and Hansen et al. (2004) reported a significant asymmetry in the dis- 
tribution of power in the WMAP data. In the reference frame defined by the north pole at 
(0, (/>) = (80°, 57°) (Galactic co-latitude and longitude), the southern hemisphere was found to 
have significantly more power on large scales {H. < 40) than the northern hemisphere. We therefore 
re-estimated the power spectra independently in these two hemispheres with the WI-FIT corrected 
maps, and found remarkably small differences with respect to the previous analyses (see Fig. 8). 
That these two different approaches to foreground subtraction yield consistent results strengthens 
the case against an entirely foreground-based explanation of the asymmetry. 

However, one flaw in this argument is that the spectral properties of the foregrounds are 
assumed to behave uniformly over the sky. We have therefore investigated this issue further by 
applying the WI-FIT procedure separately to the northern and southern hemispheres. This results 
in two sets of foreground coefficients, one for each hemisphere, that arc subsequently used to 
generate two different foreground cleaned maps. These are shown in Fig. 9 and the difference 
channel by channel is shown in Fig. 10. We see that the deviations between the maps can be quite 
large, even outside the Kp2 cut, supporting the importance of taking into account spatial variations 
of the spectral index (see the discussion in §5). However, when comparing the northern hemisphere 
power spectrum corrected by the northern coefficients with the corresponding southern hemisphere 
corrected by the southern fit, it is evident that a strong asymmetry persists, (lower plot of Fig. 8). 

As an additional check, we have applied the ILC method (Eriksen et al. 2004a) to the northern 
and southern hemisphere separately and derived consistent results. Finally, we evaluated the power 
spectrum on the northern and southern hemisphere of the map obtained by taking the frequency 
combination (2.65Ka-K)/1.65, for which synchrotron emission is strongly suppressed. Again, the 
asymmetry was very similar. If the asymmetry arises as a consequence of residual foregrounds, one 
would expect it to vary with the foregroimd subtraction procedure applied. The fact that it does 
not is a strong argument against a galactic origin for the asymmetry. 
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Table 1. WMAP band weights for the WI-FIT cleaned maps and approximate foreground 
residuals assuming the following spectral indices for the foreground components: 

/^synchrotron 2.70, /3fi-ee— free 2.15, /^thermaldust 2.20 and /^anomalousdust 2.5. The 

residuals are given relative to the expected synchrotron and anomalous dust level at 22.8 GHz, 
free-free level at 33 GHz and thermal dust level at 93.5 GHz. 



WMAP band/region 
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Ka 


Q 


V 


W 


synch. 


ff 


th. 


dust 


an. dust 


K 


1.00 


-3.75 


0.83 


3.81 


-0.90 


0.04 
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0.00 
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-1.64 
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1.00 


0.07 
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-0.06 





.46 
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W 
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Fig. 8. — The power spectra taken in the northern (lower hnes in each plot) and southern (higher 
lines in each plot) hemispheres in the frame of reference of maximum asymmetry with the north 
pole at {O^cj)) = (80°, 57°) in galactic co-latitude and longitude. Solid black lines show the results 
on the co-added ETF map obtained by the WMAP team. Green line shows the WI-FIT result for 
the Q band, blue line for the V band, red line for the W band and black crosses for the Ka band. 
In the upper plot, the results were obtained by applying the WI-FIT fitting procedure on the full 
sky (outside Kp2). In the lower plot, the fitting procedure was applied separately to the northern 
hemisphere for the northern spectra and to the southern hemisphere for the southern spectra. 
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Fig. 9. — WI-FIT cleaned WMAP maps: the maps have been produced by a noise-weighted 
combination of the separately cleaned Q, V and W bands. Upper plot shows the result of the 
fitting procedure apphed to the full sky (outside the Kp2 cut), middle/lower plot shows the same 
map with fitting applied only to the northern/southern (left/right hemisphere on the plot separated 
by a black band) hemisphere of maximum asymmetry. 




Fig. 10. — Median filtered (3 degree radius top-hat) differences between tlie noise-weighted maps 
obtained by the best fit template coefficients in the north and in the south. 
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5. Discussion and Conclusions 

We have introduced a new foreground subtraction method (WI-FIT; Wavelet based hlghreso- 
lution Fitting of Internal templates) that uses linear combinations of the foreground components, 
obtained from the microwave data itself, to fit and subtract foregrounds from the CMB dominated 
channels. For high resolution maps, the fitting procedure is performed in wavelet space since the 
pixel space approach is unfeasible for large numbers of pixels. The advantage of this method is 
that it relies neither on any assumptions nor prior knowledge about the galaxy, nor external obser- 
vations. All foreground information is obtained from the data themselves by computing differences 
between maps at different frequencies. 

We have demonstrated that the procedure works well on simulated data, based on the experi- 
mental parameters of the WMAP satellite. We then applied the method to the WMAP data and 
obtained a large-scale power spectrum in excellent agreement with that published by the WMAP 
team (see Fig. 6). We thus obtained the same results without strong constraints on the galactic 
emission based on external templates. Such a blind analysis method will be of utmost importance 
for future polarization experiments, for which reliable external templates of polarized galactic com- 
ponents may not be available. 

We further showed that the reported power spectrum asymmetry (Eriksen et al. 2004b; Hansen 
et al. 2004) remains unaltered by this approach to foreground cleaning. This holds true even when 
the foreground subtraction procedure is applied in the two opposite hemispheres separately. We 
conclude that the asymmetry is seen at the same level in all frequency bands from 23 GHz to 
94 GHz and for at least three different foreground subtraction methods. We thus consider the 
possibility of a galactic origin of the asymmetry as remote. 

In the present work, we have assumed that the frequency spectra of the foreground components 
are independent of position on the sky. This is not a very realistic assumption and also not a 
necessary one. The WI-FIT method can easily be applied to local regions on the sky individually, 
although this will likely introduce a more inhomegeneous noise structure to the data. The noisier 
areas will also need to use more large scale information in order to obtain a sufficiently high signal- 
to-noise ratio to avoid biased template coefficients, as discussed in Appendix A. The wavelet 
approach presented is well suited suited for such a scale-dependent procedure. 

Finally, we note that the the methods presented here can trivially be extended to polarization 
maps of Stoke's parameters Q and U. These issues will be investigated in a future paper. 

We acknowledge the use of the HEALPix (Gorski et al. 2005) package. FKH acknowledges 
support from a European Union Marie Curie reintegration grant. We acknowledge the use of the 
Legacy Archive for Microwave Background Data Analysis (LAMBDA). Support for LAMBDA is 
provided by the NASA Office of Space Science. 
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A. Procedures for noise bias correction and estimation of residual biases 

Because of the presence of noise in the internal templates, the template fitting procedure gives 
biased estimates of the template coefficients. This is easy to see if we take, as an illustration, fitting 
in pixel space of one noisy template in one channel. Minimizing the in equation (4), we obtain 
a best estimate 



(Al) 



Now we replace the external templates by the noisy internal templates Sj Di and write the 
internal templates in terms of a signal part and a noise part as Di = di + dni . We thus obtain 



Taking the ensemble average, we have 



(c) = c( \ (A3) 

where we have used that the CMB and the noise rij in the channel to be cleaned are both un- 
correlated with the noise Srii in the internal template. We see that when the internal template is 
noise-free, 6ni 0, the estimate is unbiased (c) c. 

We will now go to wavelet space where we will show how this bias factor can be significantly 
reduced. The bias correction procedure presented in the following can trivially be extended to real 
space. We will again start with the simplified case of one template and one channel. Then we 
extend the procedure to the realistic situation and finally we confirm that it works by showing the 
results from simulated maps. 

To simplify notation, we will now write the sum over wavelet scales as 

X*X = ^XsCss'Xs' (A4) 

SS' 

Then the best estimate can be written as 

where Xs = '^iiwiS + ^^'^ ~ Ylii'^iS ~^'^is)'^'is (^^^ equation 5) where wis, w'^g and w'^g 

are the wavelet transforms of dj, 5ni and T^" respectively. Again, the noise term is responsible 
for the bias and the first step in the bias correction procedure is to remove the noise variance term 
from Xs 

Xs^Xs- (K)^) = + 2wiwf + - (A6) 

V ' 

=Swf 
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where {{wf)^) is obtained from Monte-Carlo simulations of pure noise. The ensemble average is 
then (see equation A2) 

{6) = c{^l<£f^f±<^). (A7) 
We have noted in simulations that 

within a few percent accuracy. This suggests that we could correct for the bias if we manage to get 
the ensemble averages of the numerator and denominator in equation (A7) equal: 

The numerator and denominator can be written as 

{N) = 02 * E + * (E ^^o) (Aio) 

it i i 

(D) = {J2 ) * (E + 4((E * (E ^^o) + <(E '^^') * (E '^^')) (^n) 

i j j j i i 

Clearly, by subtracting the following terms from the denominator, the average value of the numera- 
tor and denominator will become equal and the estimate c will become unbiased provided equation 
(A8) holds. The second step of the bias correction will thus be 

D^D- 2((^ w,wf) * {J2 ^.O) - ((E '^^') * (E '^^')) (A12) 

i i i i 

The last term can be easily obtained from Monte-Carlo simulations of pure noise. The first correc- 
tion term can be obtained from simulations of noise realizations cross-correlated with the internal 
template Wi + taken from the data. Clearly how successful the bias correction is depends on 
the condition given in equation (A8) and a check of the remaining level of bias would be of high 
importance. Looking at expression (A7) it is clear that the exact expression for the bias is given 

by 

c " ^ xVx ^^^"^^ 
The problem of obtaining the bias this way is that the noise-free template Wi is unknown. We can 
however make the approximation that our noisy template is a noise-free foreground and then add 
noise realizations using Monte-Carlo simulations 



where X' = EiK+^«'^)' -2((E. O* (E. O), (x")' = ^'K+^f ) " ((E. <)*(E. O), w', = 
Wi + wf is the noisy template obtained from the data and wf are Monte-Carlo noise realizations. 
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To correct for the fact that the templates Wi and w'^ have a shghtly different structure (because of 
the noise wf), we subtract all contributions of to {N) and {D) by 

N ^ N-A{{Y,WiwfMY.'"iwf))-2{{Y,K^fMT.<'^^))-(^^ (A15) 

i i i i i i 

D ^ D-A{{Y,w^w^MY,w.w^))-^{{Y.^>f >{Y.^^^^^ (A16) 

i i i i i i 

obtained with Monte-Carlo simulations in the same manner as above. We will later show how the 
bias obtained in this way compares with the real bias in simulations. 

The bias correction procedure described above can easily be extended to the case with several 
templates and channels. Instead of making corrections to the denominator, corrections are made 
to the matrix Mff described in section §2. Using the same arguments as above, we find that the 
corresponding bias correction is given by 

Mtf ^ Mtf - Mlf - M^f - M,7, (A17) 

and the estimate of the remaining bias is obtained in the same manner as above, but again correcting 
for the additional noise term by 

Mtj ^ Mtf - Mlf - - Mlf - Mlf - Mlj - M^f - M^f - M^f - M^f, (A18) 

and 

Bf^Bf- ^(M,f + M,f + Mlf + Mff + + M/^ + M,"/), (A19) 
t 

where 



Mlf = 


( j;(^K(r) + wnT))wfit)) * EK(r) + wf{T))wf{m 

T i i 


(A20) 


Mff = 


{J2iJ2^MT) + <(r))^f (t)) * (^(^.(/) + <(/))«;f (T))) 

T i i 


(A21) 


Mff = 


{J2iJ2(Mt) + wfit))wf{T)) * {J2^Wi{T) + wf{T))wf{m 

T i i 


(A22) 


Mtf = 


T i i 


(A23) 


M^P = 




(A24) 



i i 



The matrices M^^ — M'^^ are defined similarly to M^. — M*p but with u;j — and the w'2 are drawn 
from an independent Monte-Carlo of noise realizations. Note that the used for estimating the 
bias is constructed in the same manner as the numerator in equation (A14). 

Finally these procedures were tested on WMAP-Yike simulations (described in previous sec- 
tions). In table (2) we show the size of the remaining bias on the template coefficients after applying 
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the bias correction procedure as well as estimating bias with the above method. We see that the 
bias is kept at the level of a few percent (but may be as large as 15% when considering hemispheres 
where the signal to noise is lower). Note also that the estimated bias agrees well with the real 
bias and in all cases is a conservative estimate of the bias. Note that we have not included the 
results from the W band as one of the internal templates is completely noise dominated (for the 
parameters used in the simulation) and the matrix Mff becomes unstable. 

For the WMAP data, we estimated the remaining template coefficient biases in the full sky 
fits to be less than 5% relative to the la error on the coefficients for all bands. In the hemisphere 
fits, we found the remaining bias to be in the range 5 — 20% of la. The only exception was for the 
fit to the northern hemisphere in the W-band for which the bias was estimated to be 45% of la. 
Note that the lowest wavelet scales are the scales most affected by noise. By exluding the noise 
dominated scales, the bias factor can be lowered at the cost of larger error bars. For the fit to the 
northern hemisphere in the W-band, we remade the analysis excluding the 3 first wavelet scales 
and found that the bias was reduced to a few percent. However, the results obtained with this 
map is fully consistent with the map obtained by using all scales and we have therefore chosen to 
include only the results of the latter in this paper. 
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